
rBox <- melt(df, c("survey"), c("ipl20s", "ipl16s", "ipl10s", "ipl8s"))
rBox$instr <- ifelse(grepl("20", rBox$variable), "IPL-20", 
                     ifelse(grepl("16", rBox$variable), "IPL-16",
                            ifelse(grepl("10", rBox$variable), "IPL-10", "IPL-8"
                            )
                     )
)
rBox$instr <- factor(rBox$instr, levels = c("IPL-20", "IPL-16", "IPL-10", "IPL-8"))

box10_8 <- ggplot(subset(rBox, instr == "IPL-10" | instr == "IPL-8"), aes(x = survey, y = value, fill = survey, linetype = instr)) +
  geom_boxplot(width = .4, varwidth = TRUE, size = .5, outlier.size = .3) +
  coord_cartesian(ylim = c(0, 1)) +
  scale_fill_manual(values = Col3) +
  ylab("IPL Integration Index") +
  xlab("Samples") +
  theme_minimal() +
  guides(fill = FALSE) +
  scale_x_discrete(labels = c("YouGov" = "A",
                              "Germany" = "B", 
                              "NewYork" = "C",
                              "Allies" = "D")) +
  theme(axis.title = element_text(size = 10),
        axis.text = element_text(size = 9),
        plot.title = element_text(size = 10),
        legend.position = "right")

box20_16 <- ggplot(subset(rBox, instr == "IPL-20" | instr == "IPL-16"), aes(x = survey, y = value, fill = survey, linetype = instr)) +
  geom_boxplot(width = .4, varwidth = TRUE, size = .5, outlier.size = .3) +
  coord_cartesian(ylim = c(0, 1)) +
  scale_fill_manual(values = Col3) +
  ylab("IPL Integration Index") +
  xlab("Samples") +
  theme_minimal() +
  guides(fill = FALSE) +
  scale_x_discrete(labels = c("YouGov" = "A",
                              "Germany" = "B", 
                              "NewYork" = "C",
                              "Allies" = "D")) +
  theme(axis.title = element_text(size = 10),
        axis.text = element_text(size = 9),
        plot.title = element_text(size = 10),
        legend.position = "right")

ggsave(here::here("Draft/PNAS/draft/plots", "boxpl10_8.pdf"), box10_8, height = (17.8 / 2), width = 17.8, units = "cm")
ggsave(here::here("Draft/PNAS/draft/plots", "boxpl20_16.pdf"), box20_16, height = (17.8 / 2), width = 17.8, units = "cm")

# Add t-tests

# T-test ------------------------------------------------------------------

# Create a data frame
rTests <- melt(df, c("survey"), c("ipl20s", "ipl16s", "ipl10s", "ipl8s", "ipl12s", "ipl24s"))


# Test by instrument type
rTests$measure <- ifelse(as.numeric(gsub("[[:alpha:]]+", "", rTests$variable)) > 12, "long", "short")

# Make factor for testing
rTests$group <- as.factor(paste(rTests$survey, rTests$variable, sep = ""))

# Reorder Factor labels
levels(rTests$group)
rTests$group <- factor(rTests$group,levels(rTests$group)[c(20, 19, 24, 8, 7, 12, 14, 13, 18, 2, 1, 6, 23, 22, 21, 11, 10, 9, 17, 16, 15, 5, 4, 3)])


ttpv <- pairwise.t.test(rTests$value[rTests$measure == "short"], rTests$group[rTests$measure == "short"],  p.adj = "bonf", pool.sd = FALSE)
pValuesShort <- data.frame(round(ttpv$p.value, digits = 3))
row.names(pValuesShort) <- c("A IPL-10", "A IPL-8", "B IPL-12", "B IPL-10", "B IPL-8", "C IPL-12", "C IPL-10", "C IPL-8", "D IPL-12", "D IPL-10", "D IPL-8")
names(pValuesShort) <- c("A 12", "A 10", "A 8", "B 12", "B 10", "B 8", "C 12", "C 10", "C 8", "D 12", "D 10")

ttpv <- pairwise.t.test(rTests$value[rTests$measure == "long"], rTests$group[rTests$measure == "long"],  p.adj = "bonf", pool.sd = FALSE)
pValuesLong <- data.frame(round(ttpv$p.value, digits = 3))
row.names(pValuesLong) <- c("A IPL-20", "A IPL-16", "B IPL-24", "B IPL-20", "B IPL-16", "C IPL-24", "C IPL-20", "C IPL-16", "D IPL-24", "D IPL-20", "D IPL-16")
names(pValuesLong) <- c("A 24", "A 20", "A 16", "B 24", "B 20", "B 16", "C 24", "C 20", "C 16", "D 24", "D 20")

stargazer(pValuesShort, 
          title = c("Bonferroni corrected p-values from pairwise t-tests of short measure"),
          summary = FALSE,
          rownames = TRUE,
          align = TRUE,
          digits = 3,
          font.size = "scriptsize",
          label = "pVS",
          out = c(here::here("Draft/PNAS/draft/tables", "pValuesShort.tex")))

stargazer(pValuesLong, 
          title = c("Bonferroni corrected p-values from pairwise t-tests of long measure"),
          summary = FALSE,
          rownames = TRUE,
          align = TRUE,
          digits = 3, 
          font.size = "scriptsize",
          label = "pVL",
          out = c(here::here("Draft/PNAS/draft/tables", "pValuesLong.tex")))





